home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Friedman_Rank.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-08-09  |  11.4 KB  |  639 lines
  1. /* Friedman 2-Way ANOVA by Ranks */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42.  
  43. /* Get cell reference for output range */
  44. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  45. if out_cell="" then do
  46.     'DEFPUBSCREEN("Workbench")'
  47.     exit
  48. end
  49. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  50.     'Message "Invalid cell reference"'
  51.     'DEFPUBSCREEN("Workbench")'
  52.     exit
  53. end
  54. /* Suppress Screen Redraw to Speed Things Up */
  55. 'Refresh 0'
  56.  
  57. /* Open a small output window on tcalc screen*/
  58. fo=0
  59. CR='0a'x
  60. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  61. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  62.      call writeln(6Info, DisplayMsg)
  63.     fo=1
  64. end
  65. else do
  66.     'Message "TCALC Screen not available for Progress messages"'
  67. end
  68. CALL DELAY(150)
  69.  
  70. /* Get cell references for top cell in each column */
  71. 'SelectCell' start_cell
  72. do col=start_col to end_col
  73.     'GetCursorPos'
  74.     top_cell.col=result
  75.     'Column 1'
  76. end
  77.  
  78. /* Get labels for later use on output */
  79. 'SelectCell' start_cell
  80. 'GetValue'
  81. testlabel=result
  82. testlabel=strip(testlabel)
  83. if datatype(testlabel,'n')=1 then do
  84.     labelflag=0
  85.     do x=1 to NCols
  86.     title.x="Column "||x
  87.     end
  88. end
  89. else do
  90.     labelflag=1
  91.     NRows=NRows-1
  92.     do x=1 to NCols
  93.     'GetValue'
  94.     str=result
  95.     title.x=translate(strip(str),"_"," ")
  96.     'Column 1'
  97.     end
  98. end
  99. if fo then call writech(6Info,"Progress...10 ")
  100.  
  101. /* Get data from cell range */
  102. col=start_col
  103. lav=0
  104. tot=0
  105. count.=0
  106. total.=0
  107. do x=1 to NCols
  108.     'SelectCell' top_cell.col
  109.     if labelflag=1 then 'CursorDown 1'
  110.     do y=1 to NRows
  111.         'GetValue'
  112.         valtest=result
  113.         if datatype(valtest)='NUM' then do
  114.             'GetValue'
  115.             val=result
  116.             data.x.y=val
  117.             count.x=1+count.x
  118.         end
  119.         'CursorDown 1'
  120.     end
  121.     col=col+1
  122.     val=0
  123. end
  124. if fo then call writech(6Info,"20 ")
  125.  
  126. /* Generate Sorted Rows */
  127. srow.=0
  128. Do y=1 to NRows
  129.     Do x=1 to NCols
  130.         srow.x.y=data.x.y
  131.     end
  132. call Sort()
  133. end
  134. if fo then call writech(6Info,"30 ")
  135.  
  136. /* Calculate Ranks */
  137.  
  138. Rank.=0
  139. N=NCols
  140. Do y=1 to NRows
  141.     Do x=1 to NCols
  142.     z=0
  143.         Do z=1 to N
  144.         if data.x.y=srow.z.y then Do
  145.             cnt=1
  146.             beginrank=z
  147.             avrank=z
  148.             endrank=0
  149.             Do until endrank ~=0
  150.                 z=z+1
  151.                 if (data.x.y) ~= (srow.z.y) then endrank=z-1
  152.                 avrank=avrank+z
  153.                 cnt=cnt+1
  154.                 end
  155.             val=trunc((avrank/cnt)-.5,2)
  156.             Rank.x.y=val
  157.             Leave z
  158.             end
  159.         end
  160.     end
  161. end
  162. if fo then call writech(6Info,"40 ")
  163. val=0
  164.  
  165. /* Calculate Sums etc. of Column Ranks */
  166.  
  167. sumr.=0
  168. scrsq=0
  169. Rsq.=0
  170. sr.=0
  171. A2=0
  172. B2=0
  173. Do x=1 to NCols
  174.     Do y=1 to NRows
  175.         sumr.x=(sumr.x)+(Rank.x.y)
  176.     end
  177. scrsq=scrsq+(sumr.x)**2
  178. end
  179. Do x=1 to NCols
  180.     Do y=1 to NRows
  181.         Rsq.x.y=(Rank.x.y)**2
  182.         sr.x=sr.x+Rsq.x.y
  183.     end
  184. A2=A2+sr.x
  185. end
  186. B2=scrsq/NRows
  187. if fo then call writech(6Info,"50 ")
  188.  
  189. /* Calculate Chi-Square
  190.  
  191. df=0
  192. temp=0
  193. Chi=0
  194. temp=12/((NRows*NCols)*(NCols+1))
  195. temp=temp*scrsq
  196. Chi=temp-(3*NRows)*(NCols+1)
  197. df=NCols-1*/
  198.  
  199. /* Calculate T2 Statistic */
  200.  
  201. T2=(NRows-1)*(B2-(NRows*NCols*((NCols+1)**2)/4))
  202. T2=T2/(A2-B2)
  203. k1=NCols-1
  204. k2=(NRows-1)*(NCols-1)
  205. if fo then call writech(6Info,"60 ")
  206.  
  207. /* Calculate Chi Probability
  208. P=CHIPROB(Chi,df)
  209. */
  210. /* Calculate Chi Critical
  211. Chicrit1=CHICRIT(.95,df)
  212. Chicrit2=CHICRIT(.99,DF)
  213. */
  214. /* Calculate F Probability */
  215.  
  216. AC=.0001
  217. EC=.005
  218. EC2=EC+EC
  219. PF=PROBABILITY(T2,k1,k2)
  220. PF=1-PF
  221. if fo then call writech(6Info,"70 ")
  222.  
  223. /* P=.95 OR .99 */
  224. FCRIT1=FCRITICAL(.95,k1,k2)
  225. if fo then call writech(6Info,"80 ")
  226. FCRIT2=FCRITICAL(.99,k1,k2)
  227.  
  228. if fo then do
  229.     call writeln(6Info,"100 ")
  230.     call writeln(6Info,"Writing output to window...")
  231. end
  232. /* Output */
  233. 'SelectCell' out_cell
  234. 'Column 1'
  235. 'ColumnWidth 15'
  236. 'Put "Friedman 2-Way ANOVA by Ranks"'
  237. 'CursorDown 1'
  238. 'Put "Table of Ranks"'
  239. 'CursorDown 1'
  240. Do x=1 to NCols
  241.     title=title.x
  242.     'Alignment 2'
  243.     'Put' title
  244.     'CursorDown 1'
  245.     Do y=1 to NRows
  246.         'Put' rank.x.y
  247.         'CursorDown 1'
  248.     end
  249.     'Put' sumr.x
  250. 'CursorUp' NRows+1
  251. 'Column 1'
  252. end
  253. 'SelectCell' out_cell
  254. 'CursorDown' NRows+3
  255. 'Put "Sum of Ranks:"'
  256. 'CursorDown 1'
  257. /*'Put "Chi-Square:"'
  258. 'CursorDown 1'
  259. 'Put "d.f. (v):"'
  260. 'CursorDown 1'
  261. 'Put "P(CHI<=chi):"'
  262. 'CursorDown 1'
  263. 'Put "Chi-Critical (95%):"'
  264. 'CursorDown 1'
  265. 'Put "Chi-Critical (99%):"'*/
  266. /* Put F stuff here */
  267. 'CursorDown 1'
  268. 'Put "Friedman T:"'
  269. 'CursorDown 1'
  270. 'Put "P(F<=T):"'
  271. 'CursorDown 1'
  272. 'Put "F Critical (95%):"'
  273. 'CursorDown 1'
  274. 'Put "F Critical (99%):"'
  275. 'CursorUp 3'
  276. 'Column 1'
  277. /*'Put' format(Chi,,4)
  278. 'CursorDown 1'
  279. 'Put' df
  280. 'CursorDown 1'
  281. 'Put' format(P,,6)
  282. 'CursorDown 1'
  283. 'Put' format(ChiCrit1,,4)
  284. 'CursorDown 1'
  285. 'Put' format(ChiCrit2,,4)
  286. 'CursorDown 1'*/
  287. 'Put' format(T2,,4)
  288. 'CursorDown 1'
  289. 'Put' format(PF,,6)
  290. 'CursorDown 1'
  291. 'Put' format(FCRIT1,,4)
  292. 'CursorDown 1'
  293. 'Put' format(FCRIT2,,4)
  294. 'Refresh 1'
  295. 'Refresh 2'
  296. /*'Message' "Finished"*/
  297. /*indicate the main script is finished*/
  298. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  299. result=writeln(6Info, DisplayMsg)
  300. if result~=0 then do
  301.     /*Wait 3 seconds*/
  302.     CALL DELAY(150)
  303.     /* close window*/
  304.     result=close(6Info)
  305. end
  306. 'DEFPUBSCREEN("Workbench")'
  307. exit
  308.  
  309. /* Procedures */
  310.  
  311. cellrow: procedure
  312. do
  313.     parse arg cell
  314.     do charpos=2 to length(cell)
  315.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  316.     end
  317.     return 0
  318. end
  319. Return
  320.  
  321. cellcol: procedure
  322. do
  323.     parse arg cell
  324.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  325.     cell=upper(cell)
  326.     len=length(cell)
  327.     val=0
  328. do charpos=1 to len
  329.     if datatype(substr(cell,charpos,1),n) then
  330.     do cell=reverse(substr(cell,1,charpos-1))
  331.     do x=1 to length(cell)
  332.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  333.     end
  334.     return val
  335.     end
  336.     end
  337.     return 0
  338. end
  339. Return
  340.  
  341. /* It is important to put the exposed array at the end of the next line */
  342. Sort: procedure expose y NCols srow.
  343. p=y
  344. L=(xtoy(2,int(log(NCols)/log(2))))-1
  345.     Do Until L<1
  346.     L=trunc(int(L/2))
  347.     Do J=1 to L
  348.             Do K=J+L To NCols By L
  349.             I=K
  350.             dumdat=srow.I.p
  351.             Do while I>L
  352.                 x=I-L
  353.                 If srow.x.p ~> dumdat then Leave
  354.                 srow.I.p=srow.x.p
  355.                 I=I-L
  356.             End
  357.             srow.I.p=dumdat
  358.             End
  359.         End
  360.     End
  361. Return
  362.  
  363.  
  364. syntax:
  365.      if arg(1)='FAIL' then do
  366.         'Message "Library is unavailable."'
  367.         'DEFPUBSCREEN("Workbench")'
  368.         exit
  369.         end
  370.     'DEFPUBSCREEN("Workbench")'
  371.     exit
  372.  
  373. CHIPROB: PROCEDURE
  374.  
  375.     PARSE ARG X2,K1
  376.     X3=0.5*X2
  377.     K3=0.5*K1
  378.     X=K3+1
  379.     LG=LOGGAMMA(X)
  380.     S=0
  381.     A=EXP(K3*LN(X3)-LG-X3)
  382.     IF A~=0 THEN DO
  383.         T=1
  384.         S=1
  385.         G=K3
  386.         DO WHILE (T/S) >.000001
  387.             G=G+1
  388.             T=T*X3/G
  389.             S=S+T
  390.         END
  391.     END
  392.     PO=S*A
  393. RETURN PO
  394.  
  395. CHICRIT: PROCEDURE
  396.  
  397.     PARSE ARG PO,K1
  398.     AO=.0001
  399.     E=.005
  400.     E2=E+E
  401.     PA=PO
  402.     SELECT
  403.     WHEN K1=1 THEN DO
  404.         PO=.5+PA/2
  405.         Z=NORMALPP(PO)
  406.         X1=Z*Z
  407.     END
  408.     WHEN K1=2 THEN DO
  409.         X1=1-PA
  410.         X1=-2*LN(X1)
  411.     END
  412.     OTHERWISE
  413.         Z=NORMALPP(PO)
  414.         A1=2/(9*K1)
  415.         W=1-A1+Z*SQRT(A1)
  416.         XO=K1*(W**3)
  417.         DO FOREVER
  418.             X2=XO
  419.             PO=CHIPROB(X2,K1)
  420.             F1=PO-PA
  421.             X2=XO+E
  422.             PO=CHIPROB(X2,K1)
  423.             F2=PO
  424.             X2=XO-E
  425.             PO=CHIPROB(X2,K1)
  426.             F2=F2-PO
  427.             F2=F2/E2
  428.             X1=XO-F1/F2
  429.             IF ABS(X1-XO)>AO THEN XO=X1
  430.             ELSE LEAVE
  431.         END    
  432.     END
  433.     X2=X1
  434. RETURN X2
  435.  
  436. LOGGAMMA: PROCEDURE
  437.  
  438.     ARG XA
  439.     C1=76.18009173
  440.     C2=-86.50532033
  441.     C3=24.01409822
  442.     C4=-1.231739516
  443.     C5=.001208580
  444.     C6=-.000005364
  445.     C7=2.506628275
  446.     X1=XA-1
  447.     W=X1+5.5
  448.     W=(X1+.5)*LN(W)-W
  449.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  450.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  451.     L=W+LN(C7*S)
  452. RETURN L
  453.  
  454. NORMALPP: PROCEDURE
  455.  
  456.     ARG P0
  457.     A1=2.515517
  458.     A2=.802853
  459.     A3=.010328
  460.     A4=1.432788
  461.     A5=.189269
  462.     A6=.001308
  463.     Q=.5-ABS(P0-.5)
  464.     SN=SIGN(-2*LN(Q))
  465.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  466.     ELSE W=SQRT(-2*LN(Q))
  467.     W1=A1+W*(A2+A3*W)
  468.     W2=1+W*(A4+W*(A5+A6*W))
  469.     ZZ=W-W1/W2
  470.     SELECT
  471.         WHEN (P0-.5)<0 THEN TT=-1
  472.         WHEN (P0-.5)=0 THEN TT=0
  473.         otherwise TT=1
  474.     END
  475.     ZZ=ZZ*TT
  476. RETURN ZZ
  477.  
  478. FCRITICAL: PROCEDURE EXPOSE AC EC EC2
  479.  
  480.     PARSE ARG PO,K1,K2
  481.     /* fIND NORMAL Z CORRESPONDING TO P */
  482.     P1=PO
  483.     Z=NORMALPP(PO)
  484.     IF K2<3 THEN DO
  485.         W=K2**.75
  486.         W=Z/W
  487.         Z=Z*(1.1581-W*(.2296+W*(.0042+.0027*W)))
  488.     END
  489.     /* FIND INITIAL APPROX. TO F */
  490.     A1=2/(9*K1)
  491.     A2=2/(9*K2)
  492.     W=Z*Z
  493.     W1=1+A2*(A2-W-2)
  494.     W2=A1+A2-A1*A2-1
  495.     W3=1+A1*(A1-W-2)
  496.     SN=SIGN(W2*W2-W1*W3)
  497.     IF SN=-1 THEN W3=-1*SQRT(ABS(W2*W2-W1*W3))
  498.     ELSE W3=SQRT(W2*W2-W1*W3)
  499.     FO=(W3-W2)/W1
  500.     FO=FO**3
  501.     /* MODIFIED NEWTON ITERATION TO IMPROVE VALUE OF F */
  502.     DO FOREVER
  503.         FCRIT=FO
  504.         PO=PROBABILITY(FCRIT,K1,K2)
  505.         F1=PO-P1
  506.         FCRIT=FO+EC
  507.         PO=PROBABILITY(FCRIT,K1,K2)
  508.         F2=PO
  509.         FCRIT=FO-EC
  510.         PO=PROBABILITY(FCRIT,K1,K2)
  511.         F2=F2-PO
  512.         F2=F2/EC2
  513.         F1=FO-F1/F2
  514.         IF ABS(F1-FO)>AC THEN FO=F1
  515.         ELSE LEAVE
  516.     END
  517.     FCRIT=F1
  518. RETURN FCRIT
  519.  
  520. PROBABILITY: PROCEDURE EXPOSE AC EC EC2
  521.  
  522.     PARSE ARG F,K1,K2
  523.     IF K1=1 THEN AC=.00001
  524.     H2=K2/2
  525.     H1=K1/2
  526.     H3=H1+H2
  527.     L1=0
  528.     XX=H1
  529.     LG=LOGGAMMA(XX)
  530.     L1=L1+LG
  531.     XX=H2
  532.     LG=LOGGAMMA(XX)
  533.     L1=L1+LG
  534.     XX=H3
  535.     LG=LOGGAMMA(XX)
  536.     L1=L1-LG
  537.     W1=K2/(K2+K1*F)
  538.     XX=0
  539.     IF H2<(H3*W1) THEN DO
  540.         W2=W1
  541.         W1=1-W1
  542.         W3=H2
  543.         H2=H1
  544.         H1=W3
  545.     END
  546.     ELSE DO
  547.         W2=1-W1
  548.         XX=1
  549.     END
  550.     T1=1
  551.     W3=1
  552.     P=1
  553.     I=INT(H1+W2*H3)
  554.     W4=W1/W2
  555.     /*LABELB:*/
  556.     T2=H1-W3
  557.     IF I=0 THEN W4=W1
  558.     /*LABELC:*/
  559.     DO FOREVER
  560.         T1=T1*T2*W4/(H2+W3)
  561.         P=P+T1
  562.         T2=ABS(T1)
  563.         IF (T2<=AC) & (T2<=AC*P) THEN LEAVE /*SIGNAL 'LABELD'*/
  564.         W3=W3+1
  565.         I=I-1
  566.         IF I>=0 THEN DO /*SIGNAL 'LABELB'*/
  567.             T2=H1-W3
  568.             IF I=0 THEN W4=W1
  569.             ITERATE
  570.         END
  571.         T2=H3
  572.         H3=H3+1
  573.     /*SIGNAL 'LABELC'*/
  574.     END
  575.     /*LABELD:*/
  576.     W1=EXP(H2*LN(W1)+(H1-1)*LN(W2)-L1)
  577.     P=P*W1/H2
  578.     IF XX=0 THEN PO=P
  579.     ELSE PO=1-P
  580.     AC=.001
  581. RETURN PO
  582.  
  583. Format:  procedure
  584.  
  585.      arg number, before, after
  586.      CallLine = SIGL
  587.      if ~datatype(CallLine, 'N') then CallLine = '??'
  588.  
  589.      /* Make sure we have a number as first (required) argument    */
  590.      if ~datatype(number, 'N') then do
  591.         if number = '' then
  592.            fc = 17     /* Wrong number of arguments           */
  593.         else
  594.            fc = 47     /* Arithmetic conversion error             */
  595.         signal FormatSyntaxError
  596.      end
  597.      num = number + 0
  598.      if before = '' & after = '' then
  599.         return num
  600.      else do
  601.         parse var num integer '.' fraction
  602.         if before = '' then before = length(integer)
  603.         if after = '' then after = length(fraction)
  604.         if ~datatype(before, N) | ~datatype(after, N) then
  605.            do fc = 18
  606.            signal FormatSyntaxError
  607.        end
  608.         if before < length(integer) then do
  609.            fc = 18
  610.            signal FormatSyntaxError
  611.         end
  612.         if after ~= length(fraction) then do
  613.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  614.         if integer<1&integer>-1 then integer=integer
  615.            else integer = integer + (fraction % 1)
  616.            fraction = substr(fraction, 3)
  617.         end
  618.         if fraction >= 0 then
  619.            return right(integer, before)'.'fraction
  620.         else
  621.            return right(integer, before)
  622.      end
  623.  
  624.  FormatSyntaxError:
  625.         if show('F', STDERR) then
  626.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  627.         else
  628.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  629.         'Message' mess
  630.         parse source Func .
  631.         if Func = 'FUNCTION' then do
  632.         'DEFPUBSCREEN("Workbench")'
  633.            exit "Err"
  634.         end
  635.         else do
  636.         'DEFPUBSCREEN("Workbench")'
  637.            exit 10
  638.         end
  639.